Buckling of scroll waves 
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A scroll wave in a sufficiently thin layer of an excitable medium with negative filament tension 
can be stable nevertheless due to filament rigidity. Above a certain critical thickness of the medium, 
such scroll wave will have a tendency to deform into a buckled, processing state. Experimentally 
this will be seen as meandering of the spiral wave on the surface, the amplitude of which grows 
with the thickness of the layer, until a break-up to scroll wave turbulence happens. We present a 
simplified theory for this phenomenon and illustrate it with numerical examples. 

PACS numbers: 05.45.-a, 87.23.Cc, 82.40.Ck 



Spiral waves in two-dimensions, and scroll waves in 
three-dimensions, are regimes of self-organization ob- 
served in physical, chemical and biological dissipativc 
systems, where wave propagation is supported by a 
source of energy stored in the medium [1, 2]. Due to effec- 
tive localization of the critical adjoint eigenfunctions, or 
"response functions" [3, 4], the dynamics of a spiral wave 
can be asymptotically described as that of pointwise ob- 
jects, in terms of its instant rotation centre and phase [6]. 
The third dimension endows scrolls with extra degrees of 
freedom: the filaments, around which the scroll waves ro- 
tate, can bend, and the phase of rotation may vary along 
the filaments, giving scrolls a twist [7]. The localization 
of response functions allows description of scroll waves 
as string-like objects [3, 8-12]. One manifestation of the 
extra degrees of freedom is the possibility of "scroll wave 
turbulence" due to negative tension of filaments [13]. It 
has been speculated that this scroll wave turbulence is in 
some respects similar to the hydrodynamic turbulence, 
and may provide insights into the mechanisms of cardiac 
fibrillation [3, 5, 13, 14]. 

The motivation for the present study comes from the 
analogy with hydrodynamics. At intermediate Reynolds 
numbers, laminar flow can be unstable, leading to non- 
stationary regimes which are not turbulent [2] . The pos- 
sibility of similar pre-turbulent regimes in scroll waves is 
interesting, e.g. in view of its possible relevance to cardiac 
arrhythmias. Cardiac muscle may be considered quasi- 
two-dimensional if it is very thin. Since scroll turbulence 
is essentially three-dimensional, it bears no reflection on 
behaviour of spiral waves in truly two-dimensional me- 
dia. Hence the behaviour of scrolls in layers of a given 
thickness may be effectively two-dimensional, unaffected 
by the negative tension, or truly three-dimensional, with 
full blown turbulence, or in an intermediate regime. The 
understanding of possible intermediate regimes is thus 
vitally important for interpretation of experimental data 
and for possible medical implications. 

Here we consider one such intermediate regime, which 




FIG. 1. (color online) Buckled scroll and filament, with the 
tip path on the top of the box. Barkley model with a — 1.1, 
b = 0.19, c = 0.02, D v = 0.10, box size 20 x 20 x 6.9 [15]. 



is illustrated in fig. 1. This is a snapshot of a scroll wave 
solution of an excitable reaction-diffusion model 



a t u = f(u)+DV 2 u, 



(1) 



where u,f £ M. £ , D £ M. ext , u(r,t) is the dynamic vector 
field, r £ R 3 , D is the diffusion matrix and f(u) are reac- 
tion kinetics that sustain rigidly rotating spiral waves, in 
a rectangular box r= (x, y, z) £ [0, L x ] x [0, L y ] x [0, L z ], 
with no-flux boundaries and initial conditions in the form 
of a slightly perturbed straight scroll. In boxes with L z 
below a critical height L„, the scrolls keep straight and 
rotate steadily. In large enough L z , the scroll wave tur- 
bulence develops. In a range of L z slightly above L* as 
in fig. 1, the straight scroll is unstable, and its filament, 
after an initial transient, assumes an S-like shape which 
remains constant and precesses with a constant angular 
velocity. In almost any z = const section, including the 
upper and lower surfaces, one observes spiral waves with 
a circular core, whose instant rotation centre, in turn, 



2 



rotates with an angular speed fi, which changes little 
with L z , but with a radius which is vanishingly small for 
L z }t L* and grows with L z . The resulting tip path, ob- 
served on the upper and lower surfaces, is similar to clas- 
sical two-periodic meander [16]. A similar phenomenon 
was observed in a model of heart tissue [17]. 

In this Letter, we investigate the instability which leads 
to such buckled, precessing filaments, using linear and 
non-linear theory and numerical simulations. The insta- 
bility is akin to the Euler's buckling in elasticity, where 
a straight beam deflects under a compressive stress that 
is large enough compared to the material's rigidity [18]. 

Initial insight can be obtained through linearization 
about a straight scroll wave solution U stretched along 
the z-axis, as in [24]. Small perturbations u with wave 
number k z will evolve according to 9tU = L^u, where 

Lk z = DV 2 — Dfc 2 + wade + f'(uo). (2) 

The scroll will be stable if all the eigenvalues to Lfc z have 
negative real part for all allowed wave numbers k z = 
nko = nir/L z , n G Z. Analytically, the Taylor expansion 
in k z for the critical eigenvalues A+, A_, associated to 
translational symmetry, 

A±(*„) = ±iuj -(n±ij2)k 2 z -(e 1 ±ie 2 )ki + O(kt) (3) 

relates to overlap integrals of the translational Gold- 
stone modes and response functions [3, 6, 10], see the 
Appendix [15]. With the notation of [15, 21] and tt = 
1-|V+) (W+|, wc found 

71+H2 = (W+|D|V+), (4) 

e x + ie 2 = - (W+| D(L - i^y^Ti |V+) . (5) 

Thus, a filament with negative tension 71 [3, 10, 11], can 
nevertheless be stabilized by higher order terms. We call 
e\ filament rigidity; it is an analogue of the stiffness of an 
elastic beam, and has the most important stabilizing ef- 
fect. If ei > 0, then the leading-order stability condition 
is 

k a > k* = y/— 71 /ei ^ L z < L* = iry/— ei/71. (6) 

When L z slightly exceeds L», a single unstable mode with 
spatial dependency ~ coswz/L z will grow, causing the 
filament to buckle and precess at a rate 

= 7i (71^2 - 72ei) /e\. (7) 

The amplitude at which the buckling filament will stabi- 
lize requires nonlinear analysis. Our full non-linear treat- 
ment of this phenomenon based on the time-dependent 
evolution equation for the scroll filament is rather techni- 
cal, and we defer it to another publication. Here we will 
consider simplified scroll dynamics, with the equation of 
motion of the scroll filament in the form [12] 

{9)x = (71 + 72^i?x) dlR - (ex + e 2 d a R-x) (#£)x 

+ \dlR\ 2 (b 1 + b 2 d a R*)dlR, (8) 



where R(a, t) is filament position and a is arc length. 
The coefficients bx,b 2 improve the phenomcnological rib- 
bon model proposed in [19]; they relate to the acceler- 
ated shrinking of collapsing scroll rings. Linearization of 
Eq. (8) agrees with Eqs. (6) and (7) above. A filament 
obeying Eq. (8) at L z w L* can be represented, in a in a 
frame precessing with frequency Q, by its Fourier expan- 
sion [X',Y',Z'\ = [Acos(k z),0,z} + . . . with k = tt/L z . 
Then collecting the terms ~ cos k$z gives 

A = -k 2 A [( 7l + ei kl) + k 2 A 2 q(k )] = 0, (9) 

where q(k ) = —71/2 + (3&i/4 — ex)k 2 , which describes a 
pitchfork bifurcation. By evaluating q(k*), one finds that 
the case bx > 2ei/3 yields a supercritical bifurcation, 
with stable branch 

/ 8e~x IL-L* , , 

A ^vV3^kV^' ~^ (10) 

In the opposite case, the bifurcation is subcritical. 

So, in absence of other instabilities (say two- or three- 
dimensional meander), and subject to the inequalities 
71 < 0, ex > and the limits of small I71I and small 
\L Z — L*\, we have an approximate solution (6)-(7), (10). 
The condition of negative tension, 71 < 0, is the key 
cause of the buckling instability. The condition ex > 
ensures that fourth-order arclength derivatives are suf- 
ficient to suppress high-wavenumber perturbations and 
so is important only for particular formulas but not for 
the phenomenon itself. Violation of the supercritical- 
ity condition bx > 2ei/3 does not preclude the unsta- 
ble branch from becoming stable at larger A, as will be 
seen in fig. 3(c) below. Finally, the conditions I71I <C 1, 
|Z* — Lg\ <C 1 are only required for the asymptotics; 
in reality, one would expect some finite, inter-dependent 
ranges for 71 and L z to support buckled scrolls. Hence 
we expect that buckled scrolls are fairly typical and have 
"finite chances" to be observed in some range of L z , if 
only 71 < 0. 

In our numerical simulations [15] below, the asymp- 
totics for X + (k z ) are evaluated using Eqs. (4)-(5), af- 
ter numerically obtaining the modes |V+) and (W + | 
using dxspiral [20, 21]. These asymptotics are com- 
pared to the numerical continuation of L(fc 2 )V(fc z ) = 
\ + {k z )Y{k z ) by the parameter k z [15]. 

Wc used the reaction-diffusion system (1) with 
Barkley [22] kinetics, I = 2, u = (u,v), f = (/,.g) T , 
/ = c _1 u(l — u){u — (v + b)/a), g = u — v, and 
D = diag(l, D v ). We mostly use kinetic parameters a, 6, 
c as in [23], which give negative filament tension 71 < 0, 
and consider also D v > so as to make I71 1 smaller; note 
that D v = 1 guarantees 71 = 1 > 0. 

Fig. 2(a) shows how the buckling amplitude and pre- 
cession frequency depend on the thickness of the layer, 
L z , for the same set of parameters as used to gener- 
ate fig. 1. We see that just above the bifurcation point, 
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FIG. 2. (color online) (a) Bifurcation diagram (a = 1.1, 
b = 0.19, c = 0.02, D„ = 0.1): the amplitude (upper panel) 
and precessing frequency (lower panel) of the straight and 
buckled scrolls, (b) The corresponding translational branch: 
real part (upper panel) and imaginary part (lower panel) . For 
the meandering mode, Re (A) < —0.24 [15]. 




FIG. 3. (color online) (a) Bifurcation diagram (buckling am- 
plitude) and (b) translational branch (real part), for a = 1.1, 
b = 0.19, c = 0.02, D v = 0. (c,d) Same, for a = 1.1, 
b = 0.19, c = 0.02, D v = 0.25. For the meandering modes, 
Re (A) < -0.098 and -0.32 respectively [15]. 



L z ?t L*, there is good agreement with Eq. (10). Linear 
fitting of the A 2 (L Z ) dependence for the weakest buck- 
led scrolls gives a bifurcation point L* pa 6.310, and a 
linear extrapolation of the precessing frequency from the 
same set gives Q(£*) pa 0.2789. Panel (b) shows the re- 
sults of the linear analysis, both asymptotic as given by 
Eq. (3) and obtained by numerical continuation of the 
eigenvalue problem. The latter gives the fc* pa 0.497, i.e. 
L* = ix Ik* pa 6.33, in agreement with the direct simu- 
lations shown in panel (a). The dxspiral calculations 
using Eqs. (4)-(5) give 71 = -0.353, e\ pa 2.49, result- 
ing in fc* pa 0.376. The nearly 25% difference between 
the continuation and asymptotic predictions is consis- 
tent with k z being not very small, and should decrease 
for smaller |7i|. This is indeed true, as seen below. 
The precessing frequency predicted by continuation is 
= Im(A(fc«)) - lu « 1.4188 - 1.1408 = 0.2780, in 
agreement with simulations. 

Fig. 3 illustrates variations in the buckling bifurcation 
caused by change of parameter D v . In panels (a,b) , pa- 
rameters are as in [23] and the filament tension is strongly 
negative. The dxspiral predictions are 71 pa —2.18, 
ei ~ 48.2 so the asymptotic fc* pa 0.213 is vastly dif- 
ferent from the continuation prediction fc* pa 0.890, and 
this discrepancy is clearly visible in panel (b). Yet, panel 
(a) shows that the bifurcation still takes place, and the 
critical thickness L* ps 3.60 is in agreement with the pre- 
diction of the continuation, = 7r/fc» pa 3.53. This con- 
firms that the assumption of smallness of the negative 
tension is only technical and does not preclude buckled 
scroll solutions, which still occur via a supercritical bi- 
furcation as the medium thickness varies. 



Panels (c,d) present a variation where the negative fila- 
ment tension is much smaller. Panel (d) shows much bet- 
ter agreement between the asymptotics: 71 pa —0.0362, 
ei pa 1.65, such that k* = 0.148 (L* = 21.23), whereas 
continuation gives fc* = 0.152 (L» = 20.66). However, 
the bifurcation in this case is subcritical with a hystere- 
sis, see panel (c), which shows that the assumption of 
supcrcriticality is not absolute, and that a subcritical bi- 
furcation can likewise lead to buckled scroll solutions. 

Finally, we illustrate the difference of the buckling bi- 
furcation we have described here, from the "3D meander" 
bifurcation described previously [24, 25]. On the formal 
level, the restabilized scrolls following a 3D meandering 
instability look similar: at any moment, the filament has 
a flat sinusoidal shape (given Neumann boundary con- 
ditions), and the top and bottom surfaces, as well as 
almost every z = const cross-section, show meandering 
spiral wave pictures. However, the behaviour is com- 
pletely different as L z grows, as illustrated in fig. 4. Row 
(a) shows that in the negative tension case, at sufficiently 
large L z the scroll buckles so much it breaks up and a 
scroll turbulence develops, in agreement with previous 
results. Row (b) shows that in case of 3D meander, the 
amplitude remains bounded, and even when L z is large 
enough to hold several wavelengths of the curved fila- 
ment, the restabilized "wrinkled scroll" can persist for a 
long time (compare these with "zigzag shaped filaments" 
described in [26]). Moreover, these two bifurcations occur 
in different parametric regions via different mechanisms. 
The key diffrence is, apparently, the availability or not of 
infinitely small unstable wavenumbers [15]. 
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(b) t = 300 t = 600 k z 



FIG. 4. (color online) Development of (a) autowave turbu- 
lence (a = 1.1, b = 0.19, c = 0.02, D v = 0) and (b) "wrinkled 
scroll" as restabilized solution after 3D meandering bifurca- 
tion (a = 0.66, b = 0.01, c = 0.025, D v = 0, which corre- 
sponds to the leftmost point of fig. 10(a) in [24]).Wavefronts 
are cut out by clipping planes halfway through the volume, 
to reveal the filaments. Curves on the right are real parts of 
rotational, translational and meandering eigenvalue branches 
of L fei from Eq. (2). 



To summarize, we predict that in an excitable medium 
with negative nominal filament tension 71, a sufficiently 
thin quasi-two-dimensional layer will nonetheless support 
transmural filaments which are straight and stabilized 
by filament rigidity. When the medium thickness L z 
is increased beyond a critical thickness L*, scroll waves 
may buckle and exhibit an S-shaped, processing filament. 
On the surface of the layer this will look like a classi- 
cal flower-pattern meander. If the system parameters 
yield a bifurcation of the supercritical type, a station- 
ary buckling amplitude proportional to \JL Z — will 
be reached, at which non-linear filament dynamics com- 
pensates for the negative tension 71. In the subcritical 
case, loss of stability of straight scrolls will be abrupt, 
but it still may lead to restabilized buckled scrolls. The 
knowledge about the buckling transition and its proper- 
ties is important for the planning and interpretation of 
experiments where the medium thickness is comparable 
to the typical length scale of the spiral wave. In particu- 
lar, it can be expected that stability of transmural scroll 
waves in atrial and right ventricular cardiac tissue may 
in some cases depend on filament rigidity. 
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Appendices to "Buckling of scroll waves " 
by H. Dierckx, H. Verschelde, 6. Selsil and V.N. Biktashev 



The first appendix is of theoretical nature; it clarifies 
bracket notation and offers a proof for the rigidity coef- 
ficient expression (5). The second appendix provides de- 
tails on the numerics and processing of simulation data. 
The third appendix describes extra results that might be 
of interest for some readers. Enumeration of equations 
and figures here is continued from that of the main paper, 
but the list of references is separate. 

A. Supplementary material on theory 

A.l Bracket notation for Goldstone modes and 
response functions 

We find it convenient to adopt Dirac's bra-ket notation 
from quantum mechanics and adapt them to the non- 
self adjoint problems we deal with here. Let V be a suit- 
ably chosen linear space of complex- valued m-component 
vector functions v : M 2 — > C m , the real part of which con- 
tains solutions to our reaction-diffusion system. We also 
consider its dual space, W = V*, which corresponds to 
the space of complex-valued linear functionals W[-] act- 
ing on v e V. The dual space W consists of generalized 
functions w : M 2 — > C m , which define those functionals 
via 

W[v] =JJ w H {x,y)v(x,y)dxdy = (w|v) . (11) 

So we write functions from V as ket-vectors, and func- 
tions from W as bra-vectors, assuming the scalar prod- 
uct between them when bra-vector is followed by a ket- 
vector. 

An operator A acting in V has its adjoint operator A^ 
acting in W, so that for all v € V and w 6 W, 

(AWlv) = (w|Av) , 

which is then briefly written as (w| A |v). 

In the context of spiral waves, one often linearizes the 
reaction-diffusion system (1) 

^=f(u) + DV 2 u, (12) 

around a rigidly rotating spiral wave solution U, in the 
frame that rotates with the spiral, to find 

L = DV 2 +cj a e + f'(U), 

Lt =D T V 2 -w a e + f'(U) T , 

where loq is spiral rotation rate and 9 is the polar an- 
gle. The Euclidean symmetry of the reaction-diffusion 
system (12) endows L with three critical eigenvalues 

LV (n) = A ( „ ) V(„ ) , A (n ) = inuio, n 6 {-1, 0, 1}. 



The critical eigenvectors Vm, ^(-l) which are written 
|V + ), |V_) here, are sometimes called the translational 
Goldstone modes; they can be taken in the form 

|V±) = ~\d X V±idyV). 

The spectrum of is the complex conjugate to the 
spectrum of L, and in particular 

(LtwW| = (A ( „)W(™)|, \ {n) =inu (13) 

(this is actually a nontrivial mathematical fact, see e.g. 
the discussion in [9] ) . A common choice of normalization 
is such that 

(WM|V (n) )=5- ! m,ne {-1,0,1}. (14) 

The modes (W + |, (W~| are known as 'response func- 
tions for translation'. Their belonging to W implies that 
they are effectively localized in space so that integrals 
like (11) always converge even though typical functions 
v are only bounded but not localized. Again, see [9] for 
a more detailed discussion. 

A. 2 Linearized theory 

Here we prove the result (5), which expresses the fila- 
ment rigidity coefficients ei, e2 in terms of response func- 
tions. We shall make use of a non-selfadjoint version 
of the Feynman-Hellman theorem, which states how the 
eigenvalue corresponding to a normalized eigenstate of 
a self-adjoint operator changes upon the variation of a 
real-valued parameter. For a non-selfadjoint operator A, 
if A | V) = A |V), (AtW| = (AW| and (W|V) = 1 for all 
T from a continuous interval, it follows that 

(<9 T W|V) + (W|<9 T V) = (15) 

and A = (W|AV) = (AtW|V), so 

cVA = d T (W| A | V) 

= (<9 T W| A | V) + (W| d T A | V) + (W| A |<9 T V) 
= (W|<9 T A|V). (16) 

We apply this theorem to the operator defined in Eq. (2), 
i.e. A = Lfe, = L — fc 2 D, and parametrize is by T = k z . 
We know about the continuos branches of |V(T)) and 
(W(T)|, that at k z = they reduce to the translational 
modes, |V(0)) = |V+) and (W(0)| = (W+|. 
Close to k z = 0, we expand 

A+ = iuj - jk 2 z - ek A z + 0(k 6 z ), 
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with yet unknown complex coefficients 7 = 71 + 472, 
e = e\ + iei (see Eq. (3)). In terms of the param- 
eter k\ = T, we will be looking for 7 = ~8tX+ and 
e = — ic?f,A_|_, evaluated at T = 0. From the Feynman- 
Hellman theorem (16) it follows 

&rA+ = <W| 8 T t kz \ V) = - <W| D |V) . (17) 

Evaluated at fc^ = T = 0, this recovers the well-known 
expression for the filament tension coefficient, i.e. 7 = 
7i+i72 = (W+|D|V+) [5-7]. Differentiation of Eq. (17) 
gives 

-9^A+(0) = (cVW+l D |V+) + (W+| D \d T V+) (18) 

The derivatives of the eigenfunctions can be evaluated by 
differentiating L^, |V) = A |V) with respect to T, deliv- 
ering 

(L kz - A+) \&rV) = (D d T X+) |V) , 

and similarly for (9tW|. At T = 0, we have A + = cjo, 
9tA+ = 7. We note that the linear equations for |cVV) 
and (<9tW| are solvable, because their right-hand sides 
do not have components along the null-space of the linear 
operator. Namely, it is easy to see that (D — 7) |V + ) = 
ttD |V+) and (W+| (D - 7) = (W+| Dtt, where 

7T=(1-|V + )(W + |) (19) 

is the projection operator which kills the components of 
a vector along the null space of L — iujQ. With this in 
mind, we get 

\8 T V + ) = (L — tu^fcD |V+) + C x |V+) (20) 
(cVW + | = (W+| Dtt(L - ituo)- 1 + C 2 (W+| (21) 

where C\ and C2 are arbitrary constants, which depend 
on the choices of normalizations of |V) and (W| at dif- 
ferent values of T. These choices are constrained by 
Eq. (15), which implies 

= (C 1 +C 2 )(W+|V + > 

+ (W + |D(L-iw )" 1 'r|V + ) 

+ (W+|7r(L^ r 1 D|V + ) 

= C 1 + C 2 , (22) 

because of the normalization (W + |V + ) = 1 and because 
7T |V + ) = |0) and (W+| tt = (0| by definition of n. 

Substitution of Eqs. (20) and (21) into Eq. (18), with 
account of Eq. (22), then delivers 

ei+ie 2 = --9f.A+(0) 

= (W+|D(L-^ )- 1 7rD|V + ), 

which concludes our proof of Eq. (5). Note that both 
rigidity coefficients vanish for a system with equal diffu- 
sion of variables (D = DqI), since tt |V+) = |0). 



B. Supplementary material on numerical simulations 
and data processing 

B.l Direct numerical simulations 

We used two schemes for forward evolution of the 
reaction-diffusion system, an explicit and a semi-implicit. 

Explicit scheme: Forward Euler in time with step 
At, and 7- or 19-point appoximation of the Laplacian 
with step A x in cuboid domains of size L x x L y x L z . 
We used sequential solver EZSCROLL by Barkley and 
Doyle [1], and our own sequential and MPI-parallel 
solvers. 

Semi- implicit: Operator splitting between reaction 
and diffusion substeps, with the diffusion substep by 
Brian's three-dimensional alternating-direction proce- 
dure [2, 3] which is unconditionally stable and second- 
order time accurate, implemented in our own sequential 
solver. 

The initial conditions were in the form of a straight 
scroll wave with the filament along the z coordinate, 
slightly perturbed: slightly twisted (z-dependent rota- 
tion) or slightly bended (z-dependent shift). 

The details specific for different simulation series are 
listed in Table I. The bifurcation plots in figures 2(a) 
and 3(a) each were obtained through two series of simu- 
lations: one with fixed and varied N z = L z / A x , and 
the other, to achieve a finer tuning of L Zl with fixed N z 
and varying A x . 

TABLE I. Discretization parameters for direct numerical sim- 
ulations. SI: semi-implicit; E7: explicit with 7-point Lapla- 
cian; E19: explicit with 19-point Laplacian. 



Figure 


Scheme 


A t 


A x 




L z 


1 


SI 


1/60 


1/10 


20 


6.9 


2(a) 


E7 




1/10 


16 


varied 


2(a) 


E7 


AS/12 


L z /64 


160AZ 


varied 


3(a) 


E19 


3A^/16 


1/5 


17 


varied 


3(a) 


E19 


3A^/16 


L z /19 


85A, 


varied 


3(c) 


E7 


3A^/20 


1/5 


120 


varied 


4(a,b) 


E7 


A^/12 


1/5 


40 


50 



B.2 Postprocessing of simulation data 

The results of simulations were visualized using a 
slightly modified graphical part of EZSCROLL, based on 
the Marching Cubes algorithm [1] . Figures 1 and 4 show 
snapshots of surfaces u(x,y, z,t) = u* at selected mo- 
ments of time, semi-transparent and coloured depending 
on corresponding values of v(x,y, z,t): red for smaller 
V, blue for larger v, with a smooth transition at around 
v = i>*. The tip line, which approximates the instanta- 
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neous filament, was defined as the intersection of isosur- 
faces u(x, y, z, y) = and v(x, y, z, t) = i>*, and is shown 
in green. The path of the end of the tip line at the up- 
per surface, i.e. the curve defined by u(x,y, L z ,t) = it* 
and v(x, y, L z ,t) = i>*, is drawn in grayscale, with darker 
shade corresponding to more recent position. We made 
the traditional choice for Barkley kinetics, w* = 1/2 and 
v* = a/2 — b. 

The buckling amplitude and precession were defined in 
two steps. Firstly, at a sufficiently frequent time sampling 
(t n ), say at least 30 per period, we recorded the positions 
of the tip line as X m<n = x(z m ,t n ), Y m , n = y(z m ,t n ), 
with z rn = mAj, m = 0, . . . , N z , and at each t n , approx- 
imated the tip line by 



X m>n ps A x (t„) cos(rmr/N z ), 
y m , n ~ A v (t n ) cos(mir/N z ), 



(23) 



using least squares. The resulting time series for the 
buckling amplitude vector (A x (t),A y (t)) was then aver- 
aged through periods, 



(A 



T, 



n+l 



T n +1 — T n 



A XjV (t)dt, 



with the time-integral implemented using the trapezoid 
rule. The periods were defined via u records at a selected 
point, 



U[Xr 3 Ur : Z r , Tj ) 



d t u(x r ,y r ,z r ,Tj) > 0, 



which was typically chosen in the box corner, 
(x r ,y r ,z r ) = (0,0,0). These period-averaged data were 
then used to define the amplitude A = \ (A x ) + i (A y ) | 
and phase $ = a.rg((A y ) / (A x )) of buckling. The buck- 
ling was considered established when the graph of A{t) 
showed saturation, subject to residual numerical noise. 
The value of A(t) average over a sufficiently long "es- 
tablished interval" of time was then used for graphs 
in figures 2(a) (top) and 3(a). The buckling phase was 
made "continuous" , so the difference between consecu- 
tive readings of $ does not exceed 7T, by transformation 
<f>(t) i y + 2irN t with appropriately chosen N t e Z. 
The resulting normalized dependence was approximated 
in the same established interval using least squares by a 
linear function of t, the slope of which gave the estimate 
of precession frequency fl, used for fig. 2(a) (bottom). 

For fig. 3(c), the buckling was so strong that the fil- 
ament shape was not approximated well by Eqs. (23). 
There we took instead A x (t n ) = | (Xjv^n — ^o,n)j 
A y {t n ) = | {Ypf^jTi ~ Yo.n) as the raw data. 



and (W + | were performed by dxspiral suite [8] based 
on the method described in [9], which depends on three 
discretization parameters: the disk radius /0 max , radial 
resolution N p and angular resolution Ng. For the eigen- 
value problems, we used straight shift-invert Arnoldi it- 
erations without Cayley transform, and a Krylov dimen- 
sionality of 10. The list of computed quantities used in 
previous publications [9, 10] had to be extended to com- 
pute ei + ie2, which involved the quasi- inversion process 
| a) i y |b), where 

|b)=L'7r|a) 

with L' being the inverse (L — iwo) -1 in the subspace 
orthogonal to (W+|. Recall that tt is the projection op- 
erator to that subspace, given by Eq. (19). Although the 
exact L' is not defined in the whole space, its numeri- 
cal implementation is defined, albeit extremely ill-posed. 
(For, the more accurate is the solution of the eigenvalue 
problems for |V+) and (W + |, the higher is the condition 
number of L/). Therefore, this computation presented 
some challenge. 

To achieve satisfactory results, we applied the projec- 
tion operator before and after the inverse, each several 
times, 



7T Li 7T a) 



(24) 



where k and m were integers taken as large as to ensure 
that further applications of n did not change the results 
any more at a given floating point precision (we used 
8-byte arithmetics). Obviously, the exact tt and L' com- 
mute, so mutliple applications of 7r would not change the 
result "in the ideal world" , and in the real computations 
they minimized the impact of the round-off errors and 
the magnifying effect of the ill-posed L'. 

Apart from straight application of the inverse L' 
through LU decomposition, we also tried iterative ap- 
plication of the same, a version of GMRES method and 
Tikhonov rcgularization. 

The quality of the quasi-inverse was assessed by nor- 
malized residual 

"|a)-L|b)||/|||a)||, 

where the norm is in I2. We found that with multiple 
application of 7r, the simplest method gives a satisfactory 
quality (normalized residual of the order of 10~ 2 or less) 
which is not easily improved with the other, more time- 
consuming methods. 



B.4 Continuation of the eigenvalue problem 



B.3 Numerical evaluation of the rigidity coefficients 

Computations of spiral wave solutions U, their angu- 
lar velocity loq, and their translational eigenmodes |V+) 



Our method is similar to that used in [4], up to the 
choice of the eigenvalue solver. Solving the eigenvalue 
problem 

t kz |V(fc 2 )) = X(k z ) |V(fc,)) , 
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by continuation in parameter k z , starting from a known 
initial value A(0), was done at the same discretization as 
the unperturbed spiral wave solution U and the eigen- 
modes |V + ) and (W + | of the asymptotic theory. We 
used the following discretization parameters in calculat- 
ing the eigenvalue branches: p max = 25, Ng = 1000, 
N p = 64. The problem is fully resolved at these parame- 
ters, in the sense that further increase of either of them 
docs not visibly change the graphs. 

The rotational branch Ao(fc z ) was obtained by contin- 
uation of the eigenvalue Ao(0) = 0. The translational 
branch A+(fc z ) was continued from A+(0) = iujq where 
Wo was the angular velocity of the unperturbed spiral as 
found by dxspiral. Finding the starting point for the 
for the meandering branch A m (/c z ) was more complicated. 



We have used EZRide [11] to obtain a steady spiral wave 
solution starting from cross-field initial conditions. The 
"quotient data" , representing relative velocity and angu- 
lar velocity of the tip, were approaching their equilibria 
in an oscillatory manner. After manually eliminating an 
initial transient period, these data were approximated 
by a dependency of the form Re (Ae xt ), A, A <G C, using 
Gnuplot implementation of the Marquard-Lcvenberg al- 
gorithm. Thus found A was used as an initial guess for 
the dxspiral calculations at k z = and then for contin- 
uation in k z to obtain the meandering branch. 
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C. Supplementary results 

Fig. 5 expands on the comparison of the negative ten- 
sion case leading to buckled scroll or scroll wave turbu- 
lence on one side, and the "3D meandering" instability 
leading to "wrinkled" scrolls on the other side. Here 
we have added an intermediate case (lower row), which 
shows an instability of a translational, rather than mean- 
dering, branch, however the instability is in an interval 
of k z separated from 0. The resulting phenomenology 
is the same as with 3D meandering: a seemingly stable 
wrinkled scroll is observed. Hence it appears that for the 
stability of wrinkled scrolls it is essential that the range of 
unstable wavenumbers is separated from 0, rather than 
exactly which branch shows the instability. A rigorous 
nonlinear analysis for the two cases (b) and (c) would 
have to take into account that there are several unstable 
wavenumbers in each case. 

Fig. 6 illustrates the linearization spectra for all the pa- 
rameter sets considered in the paper, shown in the same 
ranges for comparison. It is evident that spectra (a)-(d) 
show an instability of the translational mode, and spec- 
trum (c) shows an instability of the meandering mode, 
and this is not complicated by any "hybridization" de- 
scribed in [4]. The only evident hybridization is of the 
rotational mode, appearing as fracture points of the cor- 



responding Re (A) curves on panels (b) and (c). Note 
that for rotational branch Im(A) = 0. 
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